# # Quality checking of home locations and destroyed status
# # Commented because uses NearMAP API to get pre- and post-fire images of homes, so only runnable on Judd's server
# 
# library(dplyr)
# library(data.table)
# library(sf)
# library(cowplot)
# library(ggplot2)
# library(png)   #for nearmap api
# library(curl)  #for nearmap api
# library(viridis)
# 
# ## Preliminaries
# rm(list=ls())
# 
# source('code/globals.R')
# 
# 
# #------------------------------------------------------------------#
# #--- Helper functions to carry out steps of the figure making -----#
# #------------------------------------------------------------------#
# 
# ############## GET MAP TILE NUMBERS #################################################
# #### Calculate Google Map Tile x and y values for a given lat,long point and zoom
# #https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Lon..2Flat._to_tile_numbers_2
# 
# calc_tile_num<-function(lon_deg,lat_deg,zoom) {
#   ## Good notes: https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Lon..2Flat._to_tile_numbers_2
#   n<-2^zoom
#   x<-n * ((lon_deg + 180)/360)
#   y_rad<-lat_deg * (pi/180)   #https://stackoverflow.com/questions/10140029/convert-latitude-longitude-in-degree-radians
#   y <- n * (1 - (log(tan(y_rad) + (1/cos(y_rad))) / pi)) / 2    #note that 1/cos is the secant function (secant function not available in base R) 
#   xtile<-floor(x)
#   ytile<-floor(y)
#   ##-- keep the remainders -- useful for locating the point on the tile
#   xremainder<-x %% xtile
#   yremainder<-y %% ytile
#   return(c(xtile,ytile,xremainder,yremainder))
# }
# 
# ####--- Find the corner points (in lat/lng) for the tile, to georeference the raster graphic.
# ## This is the xmin/ymin of the tile, plus the xmin/ymin of the next (+1) tile (to the right and down)
# find_tile_corners<-function(tile_x,tile_y,zoom) {
#   n<-2^zoom
#   #Tile origin
#   x_this_tile = tile_x / n * 360.0 - 180.0
#   y_this_tile_rad = atan(sinh(pi * (1 - 2 * tile_y / n)))
#   y_this_tile = y_this_tile_rad * 180.0 / pi
#   rm(y_this_tile_rad)
#   #Origin for over 1 tile, down 1 tile (+1 tile in each direction)
#   x_next_tile = (tile_x+1) / n * 360.0 - 180.0
#   y_next_tile_rad = atan(sinh(pi * (1 - 2 * (tile_y+1) / n)))
#   y_next_tile = y_next_tile_rad * 180.0 / pi
#   return(c(x_this_tile,x_next_tile,y_next_tile,y_this_tile))
# }
# 
# ############ Form the API Call###########
# ## NOTE: Optional argument here to specify dates -- the default argument is Jan 1 2000 to Jan 1 2021; NM will return most recent image in that window.
# ## The argument tile should be a length-two vector with x and y tile numbers
# nm_api_call<-function(tile,zoom,apikey,
#                       time_start='2000-01-01',time_end='2021-01-01',mosaic='latest') {
#   #https://docs.nearmap.com/display/ND/Tile+API
#   call<-paste0('https://api.nearmap.com/tiles/v3/Vert/',
#                zoom,'/',tile[1],'/',tile[2],'.png?apikey=',apikey,
#                '&since=',time_start,'&until=',time_end,'&mosaic=',mosaic,'&tertiary=none')
#   return(call)
# }
# 
# gettile<-function(tile_api_call,saveloc) {
#   curl::curl_download(tile_api_call,saveloc)  # can save it somewhere better than 'dest' in the working folder
#   a<-png::readPNG(saveloc)
#   return(a)
# }
# 
# ############# Make the map figure in ggplot() ############
# # Use the png as a base_map with annotation_raster
# ## Helpful example: https://yutani.rbind.io/post/2018-06-09-plot-osm-tiles/
# ## Note for some reason, can't call ggplot() with no data -- have to put something in initial ggplot call.
# 
# make_nearmap_figure<-function(baseimage,
#                               parcelsandthelike,
#                               tiles_corners,
#                               rasterimages,
#                               testdata,
#                               bldgpolygons) {
#   
#   ##-- find the plot limits based on the tiles
#   min_xmin<-min(tiles_corners[,'xmin'])
#   max_xmax<-max(tiles_corners[,'xmax'])
#   min_ymin<-min(tiles_corners[,'ymin'])
#   max_ymax<-max(tiles_corners[,'ymax'])
#   
#   ##-- Initialize the map and choose pre-, post-imagery.
#   b<-ggplot(data=testdata)
#   ##-- add each of the 4 tiles as a georeferenced raster image using corner points
#   if (baseimage==TRUE) {
#     for (i in 1:4) {
#       b<-b+annotation_raster(rasterimages[[i]],tiles_corners[i,1],tiles_corners[i,2],
#                              tiles_corners[i,3],tiles_corners[i,4])
#     }
#   }
#   ##-- add other elements to map
#   if (parcelsandthelike==TRUE) {
#     b <- b + #geom_sf(data=testdata,color='gray',size=0.8,fill=NA) +
#       geom_sf(data=bldgpolygons,color='blue',alpha=0.1,size=0.4) +
#       #geom_sf(data=ESRI_Map_Points,aes(color=as.factor(highlightcolor)),size=1,alpha=0.5,shape=3) + #plus sign
#       geom_sf(data=testdata,aes(color=destroyed),size=1.5,alpha=0.5) + 
#       scale_color_manual(values=c("FALSE"="blue", "TRUE"="red"))
#   }
#   ##-- Add last options to map
#   b <- b + theme_map() + theme(legend.position = "none") +
#     coord_sf(xlim = c(min_xmin,max_xmax), ylim = c(min_ymin,max_ymax), expand = FALSE, crs = 4326)
#   return(b)
# }
# 
# 
# 
# #------------------------------------------------------------------#
# #--- Master function to run everything and save the figure --------#
# #------------------------------------------------------------------#
# ## Keep the nearmap option set to FALSE until ready to start downloading tiles - remember API download limits. 
# # Prototype without the base maps to find a place that looks good for checking, then turn nearmap option to TRUE.
# 
# qc_map<-function(importpcid,incid,nearmap=FALSE,zoom=18,footprintdata, lotsize=NA) {
#   
#   # importpcid<-homesinthisincident$ImportParcelID[k]
#   # incid=homesinthisincident$incidentid[k]
#   # nearmap=TRUE
#   # zoom=18
#   # lotsize=homesinthisincident$LotSizeAcres[k]
#   # footprintdata <- footprints
# 
#   ptm<-proc.time()
#   print("Setting up data and selecting building footprint polygons")
#   
#   examplehome<-merged_zillow_othervars %>% 
#     dplyr::filter(ImportParcelID==importpcid & incidentid==incid) %>%  
#     dplyr::select(lon,lat,FIPS,ImportParcelID,incidentid,State) 
#   
#   ## Get the bounding box around the example home.
#   #- because st_buffer doesn't work on lat/long data, I 
#   #- project the home location to NAD Albers, calculate the buffer,
#   #- then go back to WGS84 and find the bbox in WGS 84 coords.
#   zoombox <- st_as_sf(examplehome,coords = c("lon", "lat"),crs = 4326) %>% 
#     st_transform(crs=3310) %>% 
#     st_buffer(250) %>%   #This projection has units of meters, so this buffer should be measured in meters
#     #st_transform(crs=3857) %>% 
#     st_transform(crs=4326) %>% 
#     st_bbox()
#   
#   
#   #testdata<-st_crop(roofpoints,zoombox) %>% 
#   testdata <- merged_zillow_othervars %>% 
#     st_as_sf(coords=c("lon", "lat"), crs=4326) %>%
#     st_crop(zoombox) %>%
#     dplyr::filter(incidentid==incid) %>% ## rare cases of homes inside zip codes affected by 2 WF incidents.
#     st_transform(3857)
#   
#   sf::sf_use_s2(FALSE)
#   bldgpolygons<-st_join(footprintdata,
#                         st_transform(testdata,crs=st_crs(footprintdata)),
#                         left=FALSE) %>%
#     st_transform(crs=4326) %>%
#     #s2::s2_rebuild() %>%
#     st_crop(zoombox) %>%
#     #st_join(zoombox %>% st_as_sfc()%>%st_as_sf() %>% s2::s2_rebuild()) %>%
#     st_transform(crs=3857)
#   sf::sf_use_s2(TRUE)
#   
#   # testdata<-testdata[order(sample(nrow(testdata))),] #randomly order the rows, then color the first 10 rows
#   # testdata$accentcolor<-seq(1:nrow(testdata))
#   # if (nrow(testdata)>=11) {
#   #   testdata$accentcolor[11:nrow(testdata)]<-11
#   # }
# 
# 
#   q<-ggplot() + 
#     geom_sf(data=testdata,aes(color=destroyed),fill=NA) +
#     geom_sf(data=bldgpolygons,fill='pink',alpha=0.5) +
#     scale_color_manual(values=c("FALSE"="blue", "TRUE"="red")) + theme_map() + theme(legend.position = "none")
#   
#   lotsize_string <- ""
#   if(!is.na(lotsize)){
#     if(lotsize>1){
#       lotsize_string <- "largelot"
#     }
#     else{
#       lotsize_string <- "smalllot"
#     }
#   }
#   
# 
#   
#   mapfilename<-paste0('qcbasic_',examplehome$FIPS,"_",importpcid,"_",lotsize_string,".png")
#   save_plot(file.path(RES,"qc",mapfilename), plot = q)
#   
#   if (nearmap==TRUE) {
#     print("Getting tile information and pinging NearMap API")
#     
#     #### Get info for the tile that contains the point of interest
#     tileinfo<-calc_tile_num(examplehome[[1]],examplehome[[2]],zoom)
#     tile<-tileinfo[1:2]
#     ## Get info for the 3 adjacent tiles, on whichever side is closest to point
#     if (tileinfo[3]>0.5){
#       xmove <- 1
#     } else {
#       xmove<- -1
#     }
#     if (tileinfo[4]>0.5){
#       ymove <- 1
#     } else {
#       ymove<- -1
#     }
#     tiles<-list(tile,
#                 c(tile[1]+xmove,tile[2]),
#                 c(tile[1]+xmove,tile[2]+ymove),
#                 c(tile[1],tile[2]+ymove))
#     
#     tiles_corners<-list(find_tile_corners(tiles[[1]][1],tiles[[1]][2],zoom),
#                         find_tile_corners(tiles[[2]][1],tiles[[2]][2],zoom),
#                         find_tile_corners(tiles[[3]][1],tiles[[3]][2],zoom),
#                         find_tile_corners(tiles[[4]][1],tiles[[4]][2],zoom))
#     tiles_corners<-data.frame(xmin=c(tiles_corners[[1]][1],tiles_corners[[2]][1],tiles_corners[[3]][1],tiles_corners[[4]][1]),
#                               xmax=c(tiles_corners[[1]][2],tiles_corners[[2]][2],tiles_corners[[3]][2],tiles_corners[[4]][2]),
#                               ymin=c(tiles_corners[[1]][3],tiles_corners[[2]][3],tiles_corners[[3]][3],tiles_corners[[4]][3]),
#                               ymax=c(tiles_corners[[1]][4],tiles_corners[[2]][4],tiles_corners[[3]][4],tiles_corners[[4]][4]))
#     
#     #################################################################################
#     ############## PING NEARMAP API #################################################
#     #################################################################################
#     if (file.exists("D:\\jboomhower\\Dropbox\\Administrative\\Software\\Software\\NearMap\\NearMapAPIkey.txt")) {
#       apikey<-readLines('D:\\jboomhower\\Dropbox\\Administrative\\Software\\Software\\NearMap\\NearMapAPIkey.txt',warn=FALSE)
#     } else if (file.exists( "/home/kevin/Dropbox/UCSD/Research/NearmapAPIKey.txt")) {
#       apikey<-readLines('/home/kevin/Dropbox/UCSD/Research/NearmapAPIKey.txt', warn = FALSE)
#     }
#     
#     
#     
#     ### Load fire info
#     inc<-examplehome$incidentid
#     fire_ids<-readRDS(file.path(WORKING,"intermediatefiles/fire_ids.RDS"))
#     
#     ####-- Prefire images
#     timeperiod_start<-'2000-01-01'
#     timeperiod_end<-as.character(fire_ids$StartDate[fire_ids$incidentid==inc])
#     
#     api_calls<-lapply(tiles,nm_api_call,
#                       zoom=zoom,apikey=apikey,time_start=timeperiod_start,time_end=timeperiod_end,mosaic='latest') #mosaic = latest prioritizes most recent image within specified time window
#     print(api_calls)
#     preimages <- c()
#     successful_pre  <- tryCatch({
#       for(thistile in seq(1,4)){ # Do a loop instead of lapply to allow saving each image with a unique filename
#         preimages[[thistile]] <- gettile(api_calls[[thistile]], 
#                                        saveloc=file.path(WORKING,"intermediatefiles","qc_savedNearmapImages",paste0(incid,"_",importpcid,"_",thistile,"_pre.png")) )
#       }
#       TRUE
#     },
#      error = function(cond){
#       FALSE
#      } 
#     )
#     
#    
# 
#     ####-- Postfire images
#     timeperiod_start<-as.character(fire_ids$StartDate[fire_ids$incidentid==inc])
#     #timeperiod_end<-as.character((1.5*365)+fire_ids$StartDate[fire_ids$incidentid==inc]) #limited to 1.5 years after fire (units are days)
#     timeperiod_end <- "2021-08-19"
#     api_calls<-lapply(tiles,nm_api_call,
#                       zoom=zoom,apikey=apikey,time_start=timeperiod_start,time_end=timeperiod_end,mosaic='earliest') #mosaic=earliest will preference earliest image within the specified time window (soonest after fire, in this case)
#     #postimages<-lapply(api_calls,gettile,saveloc=file.path(WORKING,"intermediatefiles","nmimage"))
#     
#     postimages <- c()
#     successful_post <- tryCatch({
#       for(thistile in seq(1,4)){ # Do a loop instead of lapply to allow saving each image with a unique filename
#         postimages[[thistile]] <- gettile(api_calls[[thistile]], 
#                                          saveloc=file.path(WORKING,"intermediatefiles","qc_savedNearmapImages",paste0(incid,"_",importpcid,"_",thistile,"_post.png")) )
#       }
#       TRUE},
#       error=function(cond){
#         FALSE
#       }
#     )
#     #################################################################################
#     ############## MAKE THE FIGURE #################################################
#     #################################################################################
#     
#     print("Making and saving figures")
#     
#     namestem<-paste0('qcnearmap_',incid,"_",importpcid,"_",lotsize_string)
#     
#     #-- "before" imagery
#     if(successful_pre){
#       pre<-make_nearmap_figure(baseimage=TRUE,
#                                parcelsandthelike=TRUE,
#                                tiles_corners,
#                                preimages,
#                                testdata,
#                                bldgpolygons)
#       save_plot(file.path(RES,"qc",paste0(namestem,"_pre.png")), plot = pre)
#       house0$hasNearmapPre[k] <<- TRUE # assign this globally
#     }
#     #-- "after" imagery
#     if(successful_post){
#       post<-make_nearmap_figure(baseimage=TRUE,
#                                 parcelsandthelike=TRUE,
#                                 tiles_corners,
#                                 postimages,
#                                 testdata,
#                                 bldgpolygons)
#       save_plot(file.path(RES,"qc",paste0(namestem,"_post.png")), plot = post)
#       house0$hasNearmapPost[k] <<- TRUE # set this globally
#       
#       successful_maps <<- successful_maps +1 # assign this globally
#       print("Post map was successful")
#     }
#     
#     #-- no imagery
#     # noimage<-make_nearmap_figure(baseimage=FALSE,
#     #                              parcelsandthelike=TRUE,
#     #                              tiles_corners,
#     #                              preimages,
#     #                              testdata,
#     #                              bldgpolygons)
#     # save_plot(file.path(RES,"qc",paste0(namestem,"_noimage.png")), plot = noimage)
#    
#     
#     
#      #-- no parcel lines etc # MAKE ALL OF THEM WITH NO PARCEL LINES
#     # nolines<-make_nearmap_figure(baseimage=TRUE,
#     #                              parcelsandthelike=FALSE,
#     #                              tiles_corners,
#     #                              postimages,
#     #                              testdata,
#     #                              bldgpolygons)
#     # save_plot(file.path(RES,"qc",paste0(namestem,"_noparcels.png")), plot = nolines)
#     # 
#   }
#   return(paste0('Run time ',(proc.time()-ptm)[[3]] /60, ' minutes'))
# }
# 
# 
# 
# 
# #### QUALITY CHECKING EXERCISE: Drawing images randomly by incident, stratifying on lot size above/below 1 acre.
# ## Question here is API rate limits?
# 
# ## Reload zillowdata so it can be in memory without running entire script again
# zillowdata <- readRDS(file.path(WORKING,"intermediatefiles", "zillowdata.RDS"))
# 
# additional_variables<-readRDS(file.path(WORKING,"intermediatefiles", "FullRegDataDuplicatedControls.RDS")) %>% 
#   dplyr::select(ImportParcelID,BuildingOrImprovementNumber,incidentid,LotSizeAcres,IncidentAffectedHomes,totdestroyed) %>% 
#   dplyr::rename(UNIT_AND_NUMBER=incidentid)
# 
# fire_years<-readRDS(file.path(WORKING,"intermediatefiles","fire_ids.RDS")) %>% # Open this file to get the years so I can qc just the additional older fires
#   mutate(year=year(as.Date(StartDate))) %>%
#   select(year, incidentid)
# 
# merged_zillow_othervars <- zillowdata %>% 
#   merge(additional_variables,by=c('ImportParcelID','BuildingOrImprovementNumber','UNIT_AND_NUMBER'),left=TRUE) %>%
#   dplyr::rename(incidentid=UNIT_AND_NUMBER)  %>%
#   merge(fire_years, by='incidentid')
# 
# set.seed(1715121075)
# 
# # Take 10 homes from each lot size stratum in each fire 
# randomsampleimages<- merged_zillow_othervars %>% 
#   #dplyr::filter(year<=2013) %>% # Uncomment this to run for just the older fires. Has to be done this way to get the same random subset as previously generated
#   dplyr::filter(destroyed==TRUE & !is.na(LotSizeAcres)) %>% 
#   dplyr::mutate(biglot=case_when(
#     LotSizeAcres > 1 & !is.na(LotSizeAcres) ~ 1,
#     LotSizeAcres <= 1 & !is.na(LotSizeAcres) ~ 0 #,
#     #is.na(LotSizeAcres) ~ NA
#     )) %>% 
#   dplyr::group_by(incidentid,biglot) %>% 
#   dplyr::slice_sample(n=10) # Pull 10 house0s from each incident by lotsize
# 
# # Pull one home from each stratum to be "house zero"
# house0<-randomsampleimages %>% 
#   #dplyr::filter(row_number()==3) %>% # Is the third house just an arbitrary decision?
#   dplyr::filter(row_number()<=5) #%>% # Keep 5 homes per incident X lotsize
#   #dplyr::filter(totdestroyed>=20)
# 
# house0$hasNearmapPre <- FALSE
# house0$hasNearmapPost <- FALSE
# 
# 
# ### For incidents in the sample, loop over states, loading MSFT data one state at at time.
# ## Saves images that can be used for quality checking.
# statelist<-unique(house0$State)  #Why nothing in WA?
# 
# incidentlist <- unique(house0$incidentid)
# 
# 
# # random incidents to start QCing with:
# #alreadydone<- c(47, 76, 26, 67, 73, 92, 25, 95, 63, 46)
# #homeswitherrors <- c()
# 
# count <- 1
# for (incident in incidentlist){ # Start over starting with Paradise
# #incident <- incidentlist[k]
#   #if(count<=3){ # Don't do more than 3 incidents for testing
#     footprints <- file.path(WORKING, "buildings", paste0(incident,"bldgs.RDS")) %>% readRDS() %>% st_transform(crs=3857)
#     #footprints <- file.path(WORKING, "buildings", "countysubsets", paste0("b08041",".RDS")) %>% readRDS() %>% st_transform(crs=3857)
#     homesinthisincident <- house0[house0$incidentid==incident,] # pull the house0s from this state (should be 2 if there are large and small lots)
#     successful_maps <- 0 # Initialize the number of successful Nearmap API calls for this incident
# 
#     
#     if(incident %in% c("Paradise2003", "Waldo Canyon_8041_2012", "Witch2007")){# Paradise2003 and Waldo Canyon_8041_2012 have invalid geometries but this doesn't seem to work
#       footprints <- st_make_valid(footprints)
#       #homesinthisincident <- st_make_valid(homesinthisincident)
#     }
#     for (k in seq(1, nrow(homesinthisincident))){
#       if(successful_maps<2){ # make this many successful maps per incident. Give up after exhausting all house0s we kept for this incident.
#         qc_map(homesinthisincident$ImportParcelID[k],
#                 homesinthisincident$incidentid[k],
#                 nearmap=TRUE,zoom=18,
#                 footprintdata=footprints,
#                 lotsize=homesinthisincident$LotSizeAcres[k])
#             
#       }  
#     }
#     count <- count+1
#   #}
#   
# }
#  house0 %>% 
#   #filter(hasNearmapPost==TRUE || hasNearmapPre==TRUE) %>%
#   select(ImportParcelID, incidentid, FIPS, biglot, hasNearmapPre, hasNearmapPost) %>% 
#   #mutate(NoNearmap = as.numeric(ImportParcelID %in% homeswitherrors)) %>%
#   
#   write.csv(file.path(RES, "qc", "QC_AutomatedOutputOlderFires.csv"))
# 
#  
#  ####################################################################################
#  # Make Table of Statistics for manually checked QC images -----
#  ####################################################################################
#  rm(list=ls())
#  
#  source('code/globals.R')
#  
#  data <- read.csv(file.path(RES, "qc", "QC_AutomatedOutput_filled3.csv")) %>% filter(TotalCount!="NA")
#  
#  fire_ids<-readRDS(file.path(WORKING,"intermediatefiles","fire_ids.RDS"))
#  
#  neighborcounts<-readRDS(file.path(WORKING,'intermediatefiles','slimneighborvars.RDS')) %>% 
#    rename(incidentid=UNIT_AND_NUMBER) %>% 
#    select(ImportParcelID,incidentid,starts_with('totwithin')) %>% 
#    distinct(ImportParcelID,incidentid,.keep_all=T)
#  
#  loc_sources<-readRDS(file.path(WORKING,"intermediatefiles/zillowdata.RDS")) %>% 
#    select(ImportParcelID,loc_source,numfootprints) %>% 
#    distinct(ImportParcelID,.keep_all=T)
#  
#   d<-merge(data,neighborcounts,by=c('ImportParcelID','incidentid'),all.x=TRUE) %>% 
#    merge(fire_ids,by='incidentid',all.x=TRUE) %>% 
#    filter(State=="California") %>% 
#    mutate(year=year(StartDate),
#           LocationErrorRate=LocationProblems / TotalCount,
#           DamageErrorRate=DamageProblems / TotalCount) %>% 
#   merge(loc_sources,by='ImportParcelID',all.x=TRUE)
# 
#  t<-d %>% 
#    filter(!is.na(totwithin200m)) %>% 
#    mutate(bin=cut(totwithin200m,breaks=c(0,10,15),include.lowest=TRUE)) %>% 
#    group_by(bin,loc_source) %>% 
#    summarize(LocationError=weighted.mean(LocationErrorRate,w=TotalCount),
#              images=n(),
#              count=sum(TotalCount),
#              DamageError=weighted.mean(DamageErrorRate,w=TotalCount,na.rm=TRUE))
#  
# 
#  ####################################################################################
#  ### Breakdown of location sources for full dataset
#  ####################################################################################
#    
#  l<-readRDS(file.path(WORKING,"intermediatefiles/zillowdata.RDS")) %>% 
#    select(ImportParcelID,loc_source,numfootprints) %>% 
#    distinct(ImportParcelID,.keep_all=T)
#  
#  maindata<-readRDS(file.path(WORKING,"intermediatefiles/maindata.RDS")) %>% 
#    select(ImportParcelID,BuildingOrImprovementNumber,incidentid,FIPS,State) %>% 
#    merge(l,by='ImportParcelID')
#  
#  maindata %>% group_by(State,loc_source) %>% summarize(n())
#  
#  maindata %>% mutate(threeplus=numfootprints>3 & !is.na(numfootprints)) %>% 
#   group_by(threeplus,loc_source) %>% summarize(n())                      
#    
#  